5.2.1 Intro

Dal momento che la variabile dipendente pe è quantitativa continua, il problema di fitting è riconducile alla classe dei modelli supervisionati applicabili a contesti di regressione.

5.2.2 Cross-Validation

Pirma di testare i modelli è necessario: * definire la metrica da minimizzare e utilizzare per confrontare il fit del modelli: RMSE * dotarsi di un metodo di train/test che consenta di evitare il problema dell’overfitting e verificare che il modello si comporti al meglio in contesti predittivi oltre che esplicativi.

Scelgo di applicare una k-cross-validation (non ripetuta) con 5 folds: il dataset sarà diviso in 5 gruppi ed il modello sarà stimato (train) su 4 gruppi (80%) e poi applicato (test) su di 1 gruppo (20%). COn questo metodo tutte le osservazioni sono utilizzate sia per stimare il modello che per controllare il fit.

fitControl <- trainControl(method = "cv", number = 5)

5.2.3 Basic linear model

Come primo step vediamo come si comporta la regressione lineare multipla: questo modello sarà usato come benchmark per valutare le performance di modelli di ML più sofisticati:

set.seed(8848)
lm1 <- train(pe~., data = data_clean, trControl=fitControl, method = 'lm')
lm1

## Linear Regression 
## 
## 9558 samples
##    4 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 7646, 7646, 7647, 7646, 7647 
## Resampling results:
## 
##   RMSE    Rsquared   MAE     
##   4.5614  0.9287105  3.630108
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

5.2.4 ML Methods

Generalized Additive Model

# Could be slow
set.seed(8848)
gam1 <- train(pe~., data = data_clean, trControl=fitControl, method = 'gam')
gam1

## Generalized Additive Model using Splines 
## 
## 9558 samples
##    4 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 7646, 7646, 7647, 7646, 7647 
## Resampling results across tuning parameters:
## 
##   select  RMSE      Rsquared   MAE     
##   FALSE   4.161785  0.9405765  3.226169
##    TRUE   4.161912  0.9405721  3.226565
## 
## Tuning parameter 'method' was held constant at a value of GCV.Cp
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were select = FALSE and method
##  = GCV.Cp.

Boosted Generalized Additive Model

set.seed(8848)
# Could be slow
gamboost1 <- train(pe~., data = data_clean, trControl=fitControl, method = 'gamboost')
gamboost1

## Boosted Generalized Additive Model 
## 
## 9558 samples
##    4 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 7646, 7646, 7647, 7646, 7647 
## Resampling results across tuning parameters:
## 
##   mstop  RMSE      Rsquared   MAE     
##    50    4.407041  0.9344974  3.480719
##   100    4.268768  0.9376009  3.358154
##   150    4.239975  0.9383568  3.327934
## 
## Tuning parameter 'prune' was held constant at a value of no
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were mstop = 150 and prune = no.

Lasso regression

set.seed(8848)
lasso1 <- train(pe~., data = data_clean, trControl=fitControl, method = 'lasso')
lasso1

## The lasso 
## 
## 9558 samples
##    4 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 7646, 7646, 7647, 7646, 7647 
## Resampling results across tuning parameters:
## 
##   fraction  RMSE       Rsquared   MAE      
##   0.1       15.149830  0.8990594  13.144539
##   0.5        8.067583  0.9079672   6.732064
##   0.9        4.661088  0.9264860   3.717352
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was fraction = 0.9.

Ridge regression

set.seed(8848)
ridge1 <- train(pe~., data = data_clean, trControl=fitControl, method = 'ridge')
ridge1

## Ridge Regression 
## 
## 9558 samples
##    4 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 7646, 7646, 7647, 7646, 7647 
## Resampling results across tuning parameters:
## 
##   lambda  RMSE      Rsquared   MAE     
##   0e+00   4.561400  0.9287105  3.630108
##   1e-04   4.561401  0.9287104  3.630130
##   1e-01   4.873008  0.9200824  3.838092
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was lambda = 0.

Polynomial Regression Faccio variare i gradi del polinomio

# Could be slow
model_poly <- as.data.frame(matrix(0, ncol = 5, nrow = 0))
names(model_poly) <- c("i","j","k","h","RMSE")
r<-1

for(i in 1:3) {
  for(j in 1:3){
    for(k in 1:3){
      for(h in 1:3){
  f_name <-paste0("pe~poly(at,",i,")+poly(v,",j,")+poly(ap,",k,")+poly(rh,",h,")")
  form <- formula(f_name)
  set.seed(8848)
  model <- train(form, data = data_clean, trControl=fitControl, method = 'lm')
  model_poly[r,"i"] <- i
  model_poly[r,"j"] <- j
  model_poly[r,"h"] <- h
  model_poly[r,"k"] <- k
  model_poly[r,"RMSE"] <- model$results$RMSE
  r <- r+1
      }
    }
  }
}

model_poly[which.min(model_poly$RMSE),]

Linear Model with Interactions

set.seed(8848)
int1 <- train(pe~at*v*ap*rh, data = data_clean, trControl=fitControl, method = 'lm')
int1

## Linear Regression 
## 
## 9558 samples
##    4 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 7646, 7646, 7647, 7646, 7647 
## Resampling results:
## 
##   RMSE      Rsquared   MAE     
##   4.272396  0.9374021  3.337471
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

Random Forest Con CV il modello non è ottimizzato: utilizzo un approccio train-test 80-20 classico.

# Could be slow
seq <- seq(1,nrow(data_clean),1)
s <- sample(seq,as.integer(nrow(data_clean)*0.8))
data_train_rf <- data_clean[s,]
data_test_rf <- data_clean[-s,]

set.seed(8848)
rf1 <- randomForest::randomForest(pe~., data = data_train_rf, ntree = 200, mtry = 3, importance = TRUE) 
data_test_rf$predict <-predict(rf1, data_test_rf)
data_test_rf$residuals <- data_test_rf$predict - data_test_rf$pe
rf_RMSE <- sqrt(mean((data_test_rf$residuals)^2))

5.2.5 Final Model

Scelgo il modello che minimizza RMSE:

model_results <- as.data.frame(rbind(
  cbind(lm1$results$RMSE, "lm1","linear model"),
  cbind(gamboost1$results$RMSE, "gamboost1", "gamboost"),
  cbind(gam1$results$RMSE, "gam1", "gam"),
  cbind(lasso1$results$RMSE, "lasso1","lasso"),
  cbind(ridge1$results$RMSE, "ridge1","ridge"),
  cbind(model_poly[which.min(model_poly$RMSE),"RMSE"], "model_poly", "Polynomial lm"),
  cbind(int1$results$RMSE, "int1", "Lm with interactions"),
  cbind(rf_RMSE, "rf1", "random forest")
))
names(model_results) <- c("RMSE", "model", "type")
model_results$RMSE <- gsub(" ", "",model_results$RMSE)
pos_best <- which.min(as.numeric(model_results$RMSE))
model_results[pos_best,]

Summary del modello finale:

get(as.character(model_results[pos_best,"model"]))

## 
## Call:
##  randomForest(formula = pe ~ ., data = data_train_rf, ntree = 200,      mtry = 3, importance = TRUE) 
##                Type of random forest: regression
##                      Number of trees: 200
## No. of variables tried at each split: 3
## 
##           Mean of squared residuals: 10.2698
##                     % Var explained: 96.46

Importanza delle variabili:

get(as.character(model_results[pos_best,"model"]))$importance

##       %IncMSE IncNodePurity
## at 374.718235    1633981.74
## v   50.353594     496026.18
## ap  10.087611      42122.85
## rh   9.093884      37690.78